#Dan Blankenberg
#Contains file contents and helper methods for HYPHY configurations
import tempfile, os

def get_filled_temp_filename(contents):
    fh = tempfile.NamedTemporaryFile('w')
    filename = fh.name
    fh.close()
    fh = open(filename, 'w')
    fh.write(contents)
    fh.close()
    return filename

#Hard Coded hyphy path, this will need to be the same across the cluster
HYPHY_PATH = "/home/universe/linux-i686/HYPHY"
HYPHY_EXECUTABLE = os.path.join(HYPHY_PATH,"HYPHY")

NJ_tree_shared_ibf = """
COUNT_GAPS_IN_FREQUENCIES = 0;
methodIndex				  = 1;

/*-----------------------------------------------------------------------------------------------------------------------------------------*/

function InferTreeTopology(verbFlag)
{
	distanceMatrix = {ds.species,ds.species};

	MESSAGE_LOGGING 		 	= 0;
	ExecuteAFile (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"chooseDistanceFormula.def");
	InitializeDistances (0);
		
	for (i = 0; i<ds.species; i=i+1)
	{
		for (j = i+1; j<ds.species; j = j+1)
		{
			distanceMatrix[i][j] = ComputeDistanceFormula (i,j);
		}
	}

	MESSAGE_LOGGING 		 	= 1;
	cladesMade 					= 1;
	

	if (ds.species == 2)
	{
		d1 = distanceMatrix[0][1]/2;
		treeNodes = {{0,1,d1__},
					 {1,1,d1__},
					 {2,0,0}};
					 
		cladesInfo = {{2,0}};
	}
	else
	{
		if (ds.species == 3)
		{
			/* generate least squares estimates here */
			
			d1 = (distanceMatrix[0][1]+distanceMatrix[0][2]-distanceMatrix[1][2])/2;
			d2 = (distanceMatrix[0][1]-distanceMatrix[0][2]+distanceMatrix[1][2])/2;
			d3 = (distanceMatrix[1][2]+distanceMatrix[0][2]-distanceMatrix[0][1])/2;
			
			treeNodes = {{0,1,d1__},
						 {1,1,d2__},
						 {2,1,d3__}
						 {3,0,0}};
						 
			cladesInfo = {{3,0}};		
		}
		else
		{	
			njm = (distanceMatrix > methodIndex)>=ds.species;
				
			treeNodes 		= {2*(ds.species+1),3};
			cladesInfo	    = {ds.species-1,2};
			
			for (i=Rows(treeNodes)-1; i>=0; i=i-1)
			{
				treeNodes[i][0] = njm[i][0];
				treeNodes[i][1] = njm[i][1];
				treeNodes[i][2] = njm[i][2];
			}

			for (i=Rows(cladesInfo)-1; i>=0; i=i-1)
			{
				cladesInfo[i][0] = njm[i][3];
				cladesInfo[i][1] = njm[i][4];
			}
			
			njm = 0;
		}
	}
	return 1.0;
}

/*-----------------------------------------------------------------------------------------------------------------------------------------*/

function TreeMatrix2TreeString (doLengths)
{
	treeString = "";
	p = 0;
	k = 0;
	m = treeNodes[0][1];
	n = treeNodes[0][0];
	treeString*(Rows(treeNodes)*25);

	while (m)
	{	
		if (m>p)
		{
			if (p)
			{
				treeString*",";
			}
			for (j=p;j<m;j=j+1)
			{
				treeString*"(";
			}
		}
		else
		{
			if (m<p)
			{
				for (j=m;j<p;j=j+1)
				{
					treeString*")";
				}
			}	
			else
			{
				treeString*",";
			}	
		}
		if (n<ds.species)
		{
			GetString (nodeName, ds, n);
			if (doLengths != 1)
			{
				treeString*nodeName;			
			}
			else
			{	
				treeString*taxonNameMap[nodeName];
			}
		}
		if (doLengths>.5)
		{
			nodeName = ":"+treeNodes[k][2];
			treeString*nodeName;
		}
		k=k+1;
		p=m;
		n=treeNodes[k][0];
		m=treeNodes[k][1];
	}

	for (j=m;j<p;j=j+1)
	{
		treeString*")";
	}
	
	treeString*0;
	return treeString;
}
"""

def get_NJ_tree (filename):
    return """
DISTANCE_PROMPTS			= 1;
ExecuteAFile ("%s");

DataSet 				ds 			 = ReadDataFile (PROMPT_FOR_FILE);
DataSetFilter 			filteredData = CreateFilter (ds,1);

/* do sequence to branch map */

taxonNameMap = {};

for (k=0; k<ds.species; k=k+1)
{
	GetString 		(thisName, ds,k);
	shortName 		= (thisName^{{"\\\\..+",""}})&&1;
	taxonNameMap[shortName] = thisName;
	SetParameter (ds,k,shortName);
}

DataSetFilter 			filteredData = CreateFilter (ds,1);
InferTreeTopology 		(0);
treeString 				= TreeMatrix2TreeString (1);

fprintf					(PROMPT_FOR_FILE, CLEAR_FILE, treeString);
fscanf					(stdin, "String", ps_file);

if (Abs(ps_file))
{
	treeString = TreeMatrix2TreeString (2);
	UseModel (USE_NO_MODEL);
	Tree givenTree = treeString;
	baseHeight 		= TipCount (givenTree)*28;
	TREE_OUTPUT_OPTIONS = {};
	TREE_OUTPUT_OPTIONS["__FONT_SIZE__"] = 14;
	baseWidth = 0;
	treeAVL				= givenTree^0;
	drawLetter			= "/drawletter {"+TREE_OUTPUT_OPTIONS["__FONT_SIZE__"]$4+" -"+TREE_OUTPUT_OPTIONS["__FONT_SIZE__"]$2+ " show} def\\n";
	for (k3 = 1; k3 < Abs(treeAVL); k3=k3+1)
	{
		nodeName = (treeAVL[k3])["Name"];
		if(Abs((treeAVL[k3])["Children"]) == 0)
		{
			mySpecs = {};
			mySpecs ["TREE_OUTPUT_BRANCH_LABEL"] = "(" + taxonNameMap[nodeName] + ") drawLetter";
			baseWidth = Max (baseWidth, (treeAVL[k3])["Depth"]);
		}
	}
	baseWidth = 40*baseWidth;
	
	fprintf (ps_file,  CLEAR_FILE, drawLetter, PSTreeString (givenTree, "STRING_SUPPLIED_LENGTHS",{{baseWidth,baseHeight}}));
}
""" % (filename)

def get_NJ_treeMF (filename):
    return  """
ExecuteAFile ("%s");

VERBOSITY_LEVEL			= -1;
fscanf          		(PROMPT_FOR_FILE, "Lines", inLines);

_linesIn        				= Columns (inLines);
isomorphicTreesBySequenceCount  = {};

/*---------------------------------------------------------*/

_currentGene   = 1; 
_currentState = 0;
geneSeqs	  = "";
geneSeqs	  * 128;

fprintf (PROMPT_FOR_FILE, CLEAR_FILE, KEEP_OPEN);
treeOutFile = LAST_FILE_PATH;

fscanf  (stdin,"String", ps_file);
if (Abs(ps_file))
{
	fprintf (ps_file, CLEAR_FILE, KEEP_OPEN);
}

for (l=0; l<_linesIn; l=l+1)
{
	if (Abs(inLines[l]) == 0)
	{
		if (_currentState == 1)
		{
			geneSeqs      * 0;
			DataSet 	  ds 		   = ReadFromString (geneSeqs);
			_processAGene (_currentGene,treeOutFile,ps_file);
			geneSeqs * 128;
			_currentGene = _currentGene + 1;
		}
	}
	else
	{
		if (_currentState == 0)
		{
			_currentState = 1;
		}
		geneSeqs * inLines[l];
		geneSeqs * "\\n";
	}
}


if (_currentState == 1)
{
	geneSeqs      * 0;
	if (Abs(geneSeqs))
	{
		DataSet 	  ds 		   = ReadFromString (geneSeqs);
		_processAGene (_currentGene,treeOutFile,ps_file);
	}
}

fprintf (treeOutFile,CLOSE_FILE);
if (Abs(ps_file))
{
	fprintf (ps_file,CLOSE_FILE);
}
/*---------------------------------------------------------*/

function _processAGene (_geneID, nwk_file, ps_file)
{
	if (ds.species == 1)
	{
		fprintf					(nwk_file, _geneID-1, "\\tNone \\tNone\\n");
		return 0;
		
	}
	
	DataSetFilter 			filteredData = CreateFilter (ds,1);

	/* do sequence to branch map */

	taxonNameMap = {};

	for (k=0; k<ds.species; k=k+1)
	{
		GetString 		(thisName, ds,k);
		shortName 		= (thisName^{{"\\\\..+",""}});
		taxonNameMap[shortName] = thisName;
		SetParameter (ds,k,shortName);
	}

	DataSetFilter 			filteredData = CreateFilter (ds,1);
	DISTANCE_PROMPTS		= (_geneID==1);

	InferTreeTopology 		(0);
	baseTree				= TreeMatrix2TreeString (0);
	UseModel 				(USE_NO_MODEL);
	
	Tree 			baseTop	= baseTree;
	
	/* standardize this top */
	
	for (k=0; k<Abs(isomorphicTreesBySequenceCount[filteredData.species]); k=k+1)
	{
		testString	  = (isomorphicTreesBySequenceCount[filteredData.species])[k];
		Tree testTree = testString;
		if (testTree == baseTop)
		{
			baseTree = testString;
			break;
		}
	}
	if (k==Abs(isomorphicTreesBySequenceCount[filteredData.species]))
	{
		if (k==0)
		{
			isomorphicTreesBySequenceCount[filteredData.species] = {};
		}
		(isomorphicTreesBySequenceCount[filteredData.species])[k] = baseTree;
	}
	
	fprintf					(nwk_file, _geneID-1, "\\t", baseTree, "\\t", TreeMatrix2TreeString (1), "\\n");
	if (Abs(ps_file))
	{
		treeString = TreeMatrix2TreeString (2);
		UseModel (USE_NO_MODEL);
		Tree givenTree = treeString;
		baseHeight 		= TipCount (givenTree)*28;
		TREE_OUTPUT_OPTIONS = {};
		TREE_OUTPUT_OPTIONS["__FONT_SIZE__"] = 14;
		baseWidth = 0;
		treeAVL				= givenTree^0;
		drawLetter			= "/drawletter {"+TREE_OUTPUT_OPTIONS["__FONT_SIZE__"]$4+" -"+TREE_OUTPUT_OPTIONS["__FONT_SIZE__"]$2+ " show} def\\n";
		for (k3 = 1; k3 < Abs(treeAVL); k3=k3+1)
		{
			nodeName = (treeAVL[k3])["Name"];
			if(Abs((treeAVL[k3])["Children"]) == 0)
			{
				mySpecs = {};
				mySpecs ["TREE_OUTPUT_BRANCH_LABEL"] = "(" + taxonNameMap[nodeName] + ") drawLetter";
				baseWidth = Max (baseWidth, (treeAVL[k3])["Depth"]);
			}
		}
		baseWidth = 40*baseWidth;
		
		fprintf (stdout, _geneID, ":", givenTree,"\\n");
		fprintf (ps_file, PSTreeString (givenTree, "STRING_SUPPLIED_LENGTHS",{{baseWidth,baseHeight}}));
	}
	return 0;
}
""" % (filename)

BranchLengthsMF = """
VERBOSITY_LEVEL			= -1;
fscanf          		(PROMPT_FOR_FILE, "Lines", inLines);

_linesIn = Columns (inLines);

/*---------------------------------------------------------*/

_currentGene   = 1; 
_currentState = 0;
geneSeqs	  = "";
geneSeqs	  * 128;

for (l=0; l<_linesIn; l=l+1)
{
	if (Abs(inLines[l]) == 0)
	{
		if (_currentState == 1)
		{
			geneSeqs      * 0;
			DataSet 	  ds 		   = ReadFromString (geneSeqs);
			_processAGene (_currentGene);
			geneSeqs * 128;
			_currentGene = _currentGene + 1;
		}
	}
	else
	{
		if (_currentState == 0)
		{
			_currentState = 1;
		}
		geneSeqs * inLines[l];
		geneSeqs * "\\n";
	}
}

if (_currentState == 1)
{
	geneSeqs      * 0;
	if (Abs(geneSeqs))
	{
		DataSet 	  ds 		   = ReadFromString (geneSeqs);
		_processAGene (_currentGene);
	}
}

fprintf (resultFile,CLOSE_FILE);

/*---------------------------------------------------------*/

function _processAGene (_geneID)
{
	DataSetFilter 			filteredData = CreateFilter (ds,1);
	if (_currentGene == 1)
	{
		SelectTemplateModel		(filteredData);
		
		SetDialogPrompt 		("Tree file");
		fscanf	        		(PROMPT_FOR_FILE, "Tree",  givenTree);
		fscanf					(stdin, "String", resultFile);
		
		/* do sequence to branch map */
		
		validNames = {};
		taxonNameMap = {};
		
		for (k=0; k<TipCount(givenTree); k=k+1)
		{
			validNames[TipName(givenTree,k)&&1] = 1;
		}
		
		for (k=0; k<BranchCount(givenTree); k=k+1)
		{
			thisName = BranchName(givenTree,k);
			taxonNameMap[thisName&&1] = thisName;
		}
		
		storeValidNames = validNames;
		fprintf 		(resultFile,CLEAR_FILE,KEEP_OPEN,"Block\\tBranch\\tLength\\tLowerBound\\tUpperBound\\n"); 
	}
	else
	{
		validNames = storeValidNames;
	}
	
	for (k=0; k<ds.species; k=k+1)
	{
		GetString (thisName, ds,k);
		shortName = (thisName^{{"\\\\..+",""}})&&1;
		if (validNames[shortName])
		{
			taxonNameMap[shortName] = thisName;
			validNames - (shortName);
			SetParameter (ds,k,shortName);
		}
		else
		{
			fprintf 		(resultFile,"ERROR:", thisName, " could not be matched to any of the leaves in tree ", givenTree,"\\n"); 
			return 0;
		}
	}
	
	/* */
	
	LikelihoodFunction lf = (filteredData,givenTree);
	
	Optimize 				(res,lf);
	
	timer = Time(0)-timer;
	
	branchNames   = BranchName   (givenTree,-1);
	branchLengths = BranchLength (givenTree,-1);
	
	
	for (k=0; k<Columns(branchNames)-1; k=k+1)
	{
		COVARIANCE_PARAMETER = "givenTree."+branchNames[k]+".t";
		COVARIANCE_PRECISION = 0.95;
		CovarianceMatrix (cmx,lf);
		if (k==0)
		{
			/* compute a scaling factor */
			ExecuteCommands ("givenTree."+branchNames[0]+".t=1");
			scaleFactor = BranchLength (givenTree,0);
			ExecuteCommands ("givenTree."+branchNames[0]+".t="+cmx[0][1]);
		}
		fprintf (resultFile,_geneID,"\\t",taxonNameMap[branchNames[k]&&1],"\\t",branchLengths[k],"\\t",scaleFactor*cmx[0][0],"\\t",scaleFactor*cmx[0][2],"\\n");
	}
	
	ttl = (branchLengths*(Transpose(branchLengths["1"])))[0];
	global treeScaler = 1;
	ReplicateConstraint ("this1.?.t:=treeScaler*this2.?.t__",givenTree,givenTree);
	COVARIANCE_PARAMETER = "treeScaler";
	COVARIANCE_PRECISION = 0.95;
	CovarianceMatrix (cmx,lf);
	fprintf (resultFile,_geneID,"\\tTotal Tree\\t",ttl,"\\t",ttl*cmx[0][0],"\\t",ttl*cmx[0][2],"\\n");
	ClearConstraints (givenTree);
	return 0;
}
"""

BranchLengths = """
DataSet 				ds 			 = ReadDataFile (PROMPT_FOR_FILE);
DataSetFilter 			filteredData = CreateFilter (ds,1);

SelectTemplateModel		(filteredData);

SetDialogPrompt 		("Tree file");
fscanf	        		(PROMPT_FOR_FILE, "Tree",  givenTree);
fscanf					(stdin, "String", resultFile);

/* do sequence to branch map */

validNames = {};
taxonNameMap = {};

for (k=0; k<TipCount(givenTree); k=k+1)
{
	validNames[TipName(givenTree,k)&&1] = 1;
}

for (k=0; k<BranchCount(givenTree); k=k+1)
{
	thisName = BranchName(givenTree,k);
	taxonNameMap[thisName&&1] = thisName;
}

for (k=0; k<ds.species; k=k+1)
{
	GetString (thisName, ds,k);
	shortName = (thisName^{{"\\\\..+",""}})&&1;
	if (validNames[shortName])
	{
		taxonNameMap[shortName] = thisName;
		validNames - (shortName);
		SetParameter (ds,k,shortName);
	}
	else
	{
		fprintf 		(resultFile,CLEAR_FILE,"ERROR:", thisName, " could not be matched to any of the leaves in tree ", givenTree); 
		return 0;
	}
}

/* */

LikelihoodFunction lf = (filteredData,givenTree);

Optimize 				(res,lf);

timer = Time(0)-timer;

branchNames   = BranchName   (givenTree,-1);
branchLengths = BranchLength (givenTree,-1);

fprintf 		(resultFile,CLEAR_FILE,KEEP_OPEN,"Branch\\tLength\\tLowerBound\\tUpperBound\\n"); 

for (k=0; k<Columns(branchNames)-1; k=k+1)
{
	COVARIANCE_PARAMETER = "givenTree."+branchNames[k]+".t";
	COVARIANCE_PRECISION = 0.95;
	CovarianceMatrix (cmx,lf);
	if (k==0)
	{
		/* compute a scaling factor */
		ExecuteCommands ("givenTree."+branchNames[0]+".t=1");
		scaleFactor = BranchLength (givenTree,0);
		ExecuteCommands ("givenTree."+branchNames[0]+".t="+cmx[0][1]);
	}
	fprintf (resultFile,taxonNameMap[branchNames[k]&&1],"\\t",branchLengths[k],"\\t",scaleFactor*cmx[0][0],"\\t",scaleFactor*cmx[0][2],"\\n");
}

ttl = (branchLengths*(Transpose(branchLengths["1"])))[0];
global treeScaler = 1;
ReplicateConstraint ("this1.?.t:=treeScaler*this2.?.t__",givenTree,givenTree);
COVARIANCE_PARAMETER = "treeScaler";
COVARIANCE_PRECISION = 0.95;
CovarianceMatrix (cmx,lf);
ClearConstraints (givenTree);
fprintf (resultFile,"Total Tree\\t",ttl,"\\t",ttl*cmx[0][0],"\\t",ttl*cmx[0][2],"\\n");
fprintf (resultFile,CLOSE_FILE);
"""

SimpleLocalFitter = """
VERBOSITY_LEVEL               = -1;
COUNT_GAPS_IN_FREQUENCIES  = 0;

/*---------------------------------------------------------*/

function returnResultHeaders (dummy)
{
    _analysisHeaders = {};
    _analysisHeaders[0] = "GENE";
    _analysisHeaders[1] = "BP";
    _analysisHeaders[2] = "S_sites";
    _analysisHeaders[3] = "NS_sites";
    _analysisHeaders[4] = "Stop_codons";
    _analysisHeaders[5] = "LogL";
    _analysisHeaders[6] = "AC";
    _analysisHeaders[7] = "AT";
    _analysisHeaders[8] = "CG";
    _analysisHeaders[9] = "CT";
    _analysisHeaders[10] = "GT";
    _analysisHeaders[11] = "Tree";

    for (_biterator = 0; _biterator < treeBranchCount; _biterator = _biterator + 1)
    {
        branchName = treeBranchNames[_biterator];
        
        _analysisHeaders [Abs(_analysisHeaders)] = "length("+branchName+")";
        _analysisHeaders [Abs(_analysisHeaders)] = "dS("+branchName+")";
        _analysisHeaders [Abs(_analysisHeaders)] = "dN("+branchName+")";
        _analysisHeaders [Abs(_analysisHeaders)] = "omega("+branchName+")";
    }

    return _analysisHeaders;
}

/*---------------------------------------------------------*/

function runAGeneFit (myID)
{
    DataSetFilter filteredData = CreateFilter (ds,3,"","",GeneticCodeExclusions);

    if (_currentGene==1)
    {
        _MG94stdinOverload = {};
        _MG94stdinOverload ["0"] = "Local";
        _MG94stdinOverload ["1"] = modelSpecString;

        ExecuteAFile             (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"TemplateModels"+DIRECTORY_SEPARATOR+"MG94custom.mdl",
                                 _MG94stdinOverload);
        
        Tree    codonTree           = treeString;
    }
    else
    {
        HarvestFrequencies               (observedFreq,filteredData,3,1,1);
        MULTIPLY_BY_FREQS             = PopulateModelMatrix ("MG94custom", observedFreq);
        vectorOfFrequencies         = BuildCodonFrequencies (observedFreq);
        Model MG94customModel         = (MG94custom,vectorOfFrequencies,0);

        Tree    codonTree            = treeString;
    }

    LikelihoodFunction lf     = (filteredData,codonTree);

    Optimize                     (res,lf);

    _snsAVL       =                 _computeSNSSites ("filteredData", _Genetic_Code, vectorOfFrequencies, 0);
    _cL           =                  ReturnVectorsOfCodonLengths (ComputeScalingStencils (0), "codonTree");


    _returnMe = {};
    _returnMe ["GENE"]              = myID;
    _returnMe ["LogL"]              = res[1][0];
    _returnMe ["BP"]                 = _snsAVL ["Sites"];
    _returnMe ["S_sites"]             = _snsAVL ["SSites"];
    _returnMe ["NS_sites"]             = _snsAVL ["NSSites"];
    _returnMe ["AC"]                 = AC;
    _returnMe ["AT"]                 = AT;
    _returnMe ["CG"]                 = CG;
    _returnMe ["CT"]                 = CT;
    _returnMe ["GT"]                 = GT;
    _returnMe ["Tree"]                 = Format(codonTree,0,1);

    for (_biterator = 0; _biterator < treeBranchCount; _biterator = _biterator + 1)
    {
        branchName = treeBranchNames[_biterator];

        _returnMe ["length("+branchName+")"]         = (_cL["Total"])[_biterator];
        _returnMe ["dS("+branchName+")"]                 = (_cL["Syn"])[_biterator]*(_returnMe ["BP"]/_returnMe ["S_sites"]);
        _returnMe ["dN("+branchName+")"]                 = (_cL["NonSyn"])[_biterator]*(_returnMe ["BP"]/_returnMe ["NS_sites"]);
        
        ExecuteCommands ("_lom = _standardizeRatio(codonTree."+treeBranchNames[_biterator]+".nonSynRate,codonTree."+treeBranchNames[_biterator]+".synRate);");
        _returnMe ["omega("+branchName+")"]                 = _lom;
    }
    
    return _returnMe;
}

"""

SimpleGlobalFitter = """
VERBOSITY_LEVEL               = -1;
COUNT_GAPS_IN_FREQUENCIES  = 0;

/*---------------------------------------------------------*/

function returnResultHeaders (dummy)
{
    _analysisHeaders = {};
    _analysisHeaders[0]  = "GENE";
    _analysisHeaders[1]  = "BP";
    _analysisHeaders[2]  = "S_sites";
    _analysisHeaders[3]  = "NS_sites";
    _analysisHeaders[4]  = "Stop_codons";
    _analysisHeaders[5]  = "LogL";
    _analysisHeaders[6]  = "omega";
    _analysisHeaders[7]  = "omega_range";
    _analysisHeaders[8]  = "AC";
    _analysisHeaders[9]  = "AT";
    _analysisHeaders[10] = "CG";
    _analysisHeaders[11] = "CT";
    _analysisHeaders[12] = "GT";
    _analysisHeaders[13] = "Tree";

    return _analysisHeaders;
}

/*---------------------------------------------------------*/

function runAGeneFit (myID)
{
    fprintf (stdout, "[SimpleGlobalFitter.bf on GENE ", myID, "]\\n");
    taxonNameMap = {};
    
    for (k=0; k<ds.species; k=k+1)
    {
        GetString         (thisName, ds,k);
        shortName         = (thisName^{{"\\\\..+",""}})&&1;
        taxonNameMap[shortName] = thisName;
        SetParameter (ds,k,shortName);
    }

    DataSetFilter filteredData = CreateFilter (ds,1);
    _nucSites                    = filteredData.sites;
    
    if (Abs(treeString))
    {
        givenTreeString = treeString;
    }
    else
    {
        if (_currentGene==1)
        {
            ExecuteAFile             (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"Utility"+DIRECTORY_SEPARATOR+"NJ.bf");
        }
        givenTreeString = InferTreeTopology (0);
        treeString        = "";
    }

    DataSetFilter filteredData = CreateFilter (ds,3,"","",GeneticCodeExclusions);

    if (_currentGene==1)
    {
        _MG94stdinOverload = {};
        _MG94stdinOverload ["0"] = "Global";
        _MG94stdinOverload ["1"] = modelSpecString;

        ExecuteAFile             (HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"TemplateModels"+DIRECTORY_SEPARATOR+"MG94custom.mdl",
                                 _MG94stdinOverload);
        
        Tree    codonTree           = givenTreeString;
    }
    else
    {
        HarvestFrequencies               (observedFreq,filteredData,3,1,1);
        MULTIPLY_BY_FREQS             = PopulateModelMatrix ("MG94custom", observedFreq);
        vectorOfFrequencies         = BuildCodonFrequencies (observedFreq);
        Model MG94customModel         = (MG94custom,vectorOfFrequencies,0);

        Tree    codonTree            = givenTreeString;
    }

    LikelihoodFunction lf     = (filteredData,codonTree);

    Optimize                     (res,lf);

    _snsAVL       =                 _computeSNSSites ("filteredData", _Genetic_Code, vectorOfFrequencies, 0);
    _cL           =                  ReturnVectorsOfCodonLengths (ComputeScalingStencils (0), "codonTree");


    _returnMe = {};
    _returnMe ["GENE"]              = myID;
    _returnMe ["LogL"]              = res[1][0];
    _returnMe ["BP"]                 = _snsAVL ["Sites"];
    _returnMe ["S_sites"]             = _snsAVL ["SSites"];
    _returnMe ["NS_sites"]             = _snsAVL ["NSSites"];
    _returnMe ["Stop_codons"]         = (_nucSites-filteredData.sites*3)$3;
    _returnMe ["AC"]                 = AC;
    _returnMe ["AT"]                 = AT;
    _returnMe ["CG"]                 = CG;
    _returnMe ["CT"]                 = CT;
    _returnMe ["GT"]                 = GT;
    _returnMe ["omega"]             = R;
    COVARIANCE_PARAMETER             = "R";
    COVARIANCE_PRECISION             = 0.95;
    CovarianceMatrix                 (cmx,lf);
    _returnMe ["omega_range"]         = ""+cmx[0]+"-"+cmx[2];
    _returnMe ["Tree"]                 = Format(codonTree,0,1);


    return _returnMe;
}
"""

FastaReader = """
fscanf			(stdin, "String", _coreAnalysis);
fscanf			(stdin, "String", _outputDriver);

ExecuteAFile 			(HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"TemplateModels"+DIRECTORY_SEPARATOR+"chooseGeneticCode.def");
ExecuteAFile 			(HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"dSdNTreeTools.ibf");
ExecuteAFile 			(HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"Utility"+DIRECTORY_SEPARATOR+"CodonTools.bf");
ExecuteAFile 			(HYPHY_BASE_DIRECTORY+"TemplateBatchFiles"+DIRECTORY_SEPARATOR+"Utility"+DIRECTORY_SEPARATOR+"GrabBag.bf");

SetDialogPrompt ("Tree file");
fscanf	        (PROMPT_FOR_FILE, "Tree",  givenTree);

treeBranchNames			   = BranchName (givenTree,-1);
treeBranchCount			   = Columns    (treeBranchNames)-1;
treeString 				   = Format (givenTree,1,1);

SetDialogPrompt ("Multiple gene FASTA file");
fscanf          (PROMPT_FOR_FILE, "Lines", inLines);
fscanf			(stdin, "String",  modelSpecString);
fscanf			(stdin, "String", _outPath);
 
ExecuteAFile	(_outputDriver);
ExecuteAFile	(_coreAnalysis);
 
/*---------------------------------------------------------*/

_linesIn     = Columns (inLines);
_currentGene   = 1; 
 _currentState = 0;
/* 0 - waiting for a non-empty line */
/* 1 - reading files */

geneSeqs       = "";
geneSeqs 	   * 0;

_prepareFileOutput (_outPath);

for (l=0; l<_linesIn; l=l+1)
{
	if (Abs(inLines[l]) == 0)
	{
		if (_currentState == 1)
		{
			geneSeqs      * 0;
			DataSet 	  ds 		   = ReadFromString (geneSeqs);
			_processAGene (ds.species == treeBranchCount,_currentGene);
			geneSeqs * 128;
			_currentGene = _currentGene + 1;
		}
	}
	else
	{
		if (_currentState == 0)
		{
			_currentState = 1;
		}
		geneSeqs * inLines[l];
		geneSeqs * "\\n";
	}
}

if (_currentState == 1)
{
	geneSeqs      * 0;
	DataSet 	  ds 		   = ReadFromString (geneSeqs);
	_processAGene (ds.species == treeBranchCount,_currentGene);
}

_finishFileOutput (0);
"""

TabWriter = """
/*---------------------------------------------------------*/
function _prepareFileOutput (_outPath)
{
	_outputFilePath = _outPath;
	
	_returnHeaders 	= returnResultHeaders(0);
	
	fprintf (_outputFilePath, CLEAR_FILE, KEEP_OPEN, _returnHeaders[0]);
	for (_biterator = 1; _biterator < Abs(_returnHeaders); _biterator = _biterator + 1)
	{
		fprintf (_outputFilePath,"\\t",_returnHeaders[_biterator]);
	}


	
	fprintf (_outputFilePath,"\\n");
	return 0;
}

/*---------------------------------------------------------*/

function _processAGene (valid, _geneID)
{
	if (valid)
	{
		returnValue = runAGeneFit (_geneID);
		fprintf (_outputFilePath, returnValue[_returnHeaders[0]]);
		for (_biterator = 1; _biterator < Abs(_returnHeaders); _biterator = _biterator + 1)
		{
			fprintf (_outputFilePath,"\\t",returnValue[_returnHeaders[_biterator]]);
		}
		fprintf (_outputFilePath, "\\n");
	}
	else
	{
		fprintf (_outputFilePath, 
				_geneID, ", Incorrect number of sequences\\n");
	}
	_currentState = 0;
	return 0;
}

/*---------------------------------------------------------*/
function _finishFileOutput (dummy)
{
	return 0;
}
"""

def get_dnds_config_filename(Fitter_filename, TabWriter_filename, genetic_code, tree_filename, input_filename, nuc_model, output_filename, FastaReader_filename ):
    contents = """
_genomeScreenOptions = {};

/* all paths are either absolute or relative 
to the DATA READER */

_genomeScreenOptions ["0"] = "%s"; 
    /* which analysis to run on each gene; */    
_genomeScreenOptions ["1"] = "%s"; 
    /* what output to produce; */    

_genomeScreenOptions ["2"] = "%s"; 
    /* genetic code */    
_genomeScreenOptions ["3"] = "%s"; 
    /* tree file */
_genomeScreenOptions ["4"] = "%s"; 
    /* alignment file */
_genomeScreenOptions ["5"] = "%s"; 
    /* nucleotide bias string; can define any of the 203 models */
_genomeScreenOptions ["6"] = "%s"; 
    /* output csv file */
    
ExecuteAFile ("%s", _genomeScreenOptions);    
""" % (Fitter_filename, TabWriter_filename, genetic_code, tree_filename, input_filename, nuc_model, output_filename, FastaReader_filename )
    return get_filled_temp_filename(contents)
    
    
def get_branch_lengths_config_filename(input_filename, nuc_model, model_options, base_freq, tree_filename, output_filename, BranchLengths_filename):
    contents = """
_genomeScreenOptions = {};

/* all paths are either absolute or relative 
to the NucDataBranchLengths.bf */

_genomeScreenOptions ["0"] = "%s"; 
	/* the file to analyze; */	
_genomeScreenOptions ["1"] = "CUSTOM"; 
	/* use an arbitrary nucleotide model */	
_genomeScreenOptions ["2"] = "%s";
	/* which model to use */ 
_genomeScreenOptions ["3"] = "%s"; 
	/* model options */
_genomeScreenOptions ["4"] = "Estimated"; 
	/* rate parameters */
_genomeScreenOptions ["5"] = "%s"; 
	/* base frequencies */
_genomeScreenOptions ["6"] = "%s"; 
	/* the tree to use; */	
_genomeScreenOptions ["7"] = "%s"; 
	/* write .csv output to; */	
	
ExecuteAFile ("%s", _genomeScreenOptions);	
""" % (input_filename, nuc_model, model_options, base_freq, tree_filename, output_filename, BranchLengths_filename)
    return get_filled_temp_filename(contents)
    
    
def get_nj_tree_config_filename(input_filename, distance_metric, output_filename1, output_filename2, NJ_tree_filename):
    contents = """
_genomeScreenOptions = {};

/* all paths are either absolute or relative 
to the BuildNJTree.bf */

_genomeScreenOptions ["0"] = "%s"; 
	/* the file to analyze; */	
_genomeScreenOptions ["1"] = "%s"; 
	/* pick which distance metric to use; TN93 is a good default */	
_genomeScreenOptions ["2"] = "%s"; 
	/* write Newick tree output to; */	
_genomeScreenOptions ["3"] = "%s"; 
	/* write a postscript tree file to this file; leave blank to not write a tree */	
	
ExecuteAFile ("%s", _genomeScreenOptions);	
""" % (input_filename, distance_metric, output_filename1, output_filename2, NJ_tree_filename)
    return get_filled_temp_filename(contents)
    
    
def get_nj_treeMF_config_filename(input_filename, output_filename1, output_filename2, distance_metric, NJ_tree_filename):
    contents = """
_genomeScreenOptions = {};

/* all paths are either absolute or relative 
to the BuildNJTreeMF.bf */

_genomeScreenOptions ["0"] = "%s"; 
	/* the multiple alignment file to analyze; */	
_genomeScreenOptions ["1"] = "%s"; 
	/* write Newick tree output to; */	
_genomeScreenOptions ["2"] = "%s"; 
	/* write a postscript tree file to this file; leave blank to not write a tree */	
_genomeScreenOptions ["3"] = "%s"; 
	/* pick which distance metric to use; TN93 is a good default */	
	
ExecuteAFile ("%s", _genomeScreenOptions);		
""" % (input_filename, output_filename1, output_filename2, distance_metric, NJ_tree_filename)
    return get_filled_temp_filename(contents)